Orbital degrees of freedom as origin of magnetoelectric coupling in magnetite 
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A microscopic understanding of magnetoelectricity, i. e. the coupling between magnetic (electric) 
properties and external electric (magnetic) fields, is a crucial milestone for future generations of 
electrically-controlled spintronic devices. Here, we focus on the first magnetoelectric known to 
mankind: magnetite. By means of a joint approach based on phenomenological Landau theory and 
density- functional simulations, we show that magnetoelectricity in charge-/orbital-ordered Fe^O^ in 
the non-centrosymmetric Cc structure is driven by the interplay between a peculiar orbital-order 
and on-site spin-orbit coupling. The excellent agreement with available experiments confirms our 
theoretical picture, pointing to magnetite as a prototype of a novel category of magnetoelectrics 
where ferroelectric polarization can be induced, tuned or switched via a magnetic field. 

PACS numbers: Valid PACS appear here 



Magnetoelectric (ME) effects - - i.e. how to con- 
trol magnetic (electric) properties via electric (mag- 
netic) fields — have been intensively investigated in re- 
cent years, often in connection to multiferroics, mate- 
rials where two or more ferroic properties (i.e. long 
range spontaneous magnetism, deformation or ferroelec- 
tricity) coexist in the same phase. [H-Q In this framework, 
it was recently found ferrimagnetic magnetite (FesC^) 
to show a sizable ferroelectric (FE) polarization at low 
temperatures .with the microscopic mechanism clarified 
theoretically: [4, 5] the FE polarization, of the order of few 
/iC/cm 2 , is induced purely by the non-centrosymmetric 
charge-order (CO), a mechanism that is rarely seen in 
other materials. Q 

The ME effect in magnetite was discovered by Rado 
and Ferrari in 1975 [9]; however, its importance has been 
neglected for long time. It is the aim of this paper to 
clarify the microscopic origin of this peculiar ME effect 
in magnetite, by means of Landau theory and density- 
functional-theory (DFT). We start from two important 
observations made in Ref. S (i) The electric polariza- 
tion P behaves as a nonlinear function of the applied 
magnetic field H, more precisely, P a oc sin 2 (with neg- 
ligible sin 4 # term), with denoting the polar angle of 
the H direction at 4.2K. (ii) The observed P (reported 
in arbitrary units) is saturated with respect to a certain 
i7 sat , i.e. its behavior being similar to the magnetization 
M. This implies that P arises in the perfectly "collinear" 
ferrimagnetic configuration, with all spins simultaneously 
following the direction of H (i.e. M / /H); this shows a 
distinct contrast with respect to many cases of single- 
phase magnetoelectrics, where either ii"-induced "non- 
collinear" spin canting or spin spiral modulation occur 
(e.g. at Cr 2 O 3 [10], LiNiP0 4 [lll and TbMn0 3 [H). 

Structural Details. — Magnetite shows a well-known 
metal-insulator transition, namely the "Verwey tran- 
sition" at Tv=120K, below which the crystal struc- 
ture changes from cubic Fd3m to a less symmetric 
structure; correspondingly, the Fe 2+ /Fe 3+ charge or- 
dering, observed on Fe-B sites in the inverse-spinel 
AB2O4 lattice, [3, [3] arises. Despite several differ- 




FIG. 1: (a) Charge/orbital ordering shown by an isosurface of 
charge density of Fe-t^ states in the upper half (l/2<z< 1) 
of the primitive Cc unit cell (for details, see ref. i). The 
red arrow shows the electric dipole connecting B12(Fe 2+ ) and 
B14(Fe 3+ ) sites, responsible for the ferroelectric P. Vectors 
(a, 6, c) denote directions of the conventional lattice vectors, 
linked to the cartesian coordinate system through x=a + 6, 
y=—a + b and z ~ c (due to a small monoclinic distortion). 
The occupied orbital state at some Fe 2+ sites is labeled as 
xy/yz/zx, used as site index hereafter, (b) Change of the 
charge density (an isosurface and a section) on Fe-B 12 site 
upon rotation of the spin direction (shown by blue arrows 
in the upper panel). To emphasize the effect, SOC term is 
enhanced by a factor of 10. 

ent space groups were proposed for the low-temperature 
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structure, [ 14j a base-centered monoclinic Cc symmetry 
(#9) was suggested as ground state both by experimen- 
tal and theoretical studies. [HI, Gil The absence of inver- 
sion symmetry in the Cc space group and the consequent 
polar CO pattern were found to be relevant for the pres- 
ence of ferroelectricity. [4] Indeed, the FE polarization is 
primarily caused by uncompensated local electric dipoles 
connecting Fe 2+ and Fe 3+ ions at B sites; as shown in 
FigUJa) and deeply discussed in Refs.[E|,[l7|, Fe 2+ at B12 
site and Fe 3+ ions at B14 site form one of the dipoles re- 
sponsible for the net P in the Cc unit cell. FigJT^a) also 
shows the orbital order (00) at Fe 2+ (d 6 ) sites. Being 
coupled to cooperative Jahn- Teller (JT) distortions, the 
partially filled minority-spin t\^ g orbital shows one of the 
three d xy /d yz /d zx charge distributions. In the primitive 
unit cell, 32 Fe-B sites are split into 16 Fe 3+ sites and 
16 Fe 2+ sites; in addition, the latter JT-active sites are 
furthermore split into three groups: 8 xy, 4 yz and 4 
zx sites. [18| The 00 pattern is such that the orbitals 
avoid to overlap, consistent with reducing the inter-site 
Coulomb repulsion and optimizing cooperative JT dis- 
tortions. Consequently, these orbitals do not fully lie in 
the xy/yz/zx planes, but are rather slightly tilted. 

Macroscopic Model. - Hereafter, we separate the 
spontaneous CO-induced, P co , and the M-direction- 
dependent, P ME (A4"). According to Ref. [TH, the mecha- 
nisms leading to P ME are generally classified as driven by 
magnetostriction (P MS ), spin-current (P sp ) and orbital 
degrees of freedom (P orb ). In our case, the collinear fer- 
rimagnetic configuration rules out both P MS and P sp , 
whereas only P orb can be induced by the spin-orbit cou- 
pling (SOC) within orbital-unquenched t\ g e? g state at 
Fe 2+ (d 6 ) sites. Summarizing, the total polarization is 
here described as P total (M) = P co + P orb (M). As 
we will discuss later, the P orb term is deeply related to 
the single-site magnetic anisotropy energy (MAE), both 
having SOC as common underlying origin. 

TABLE I: Matrices of the generators of the Cc space group 
in the representations spanned by M and P. The space group 
elements are denoted as the identity E and c-glide={crfe|00|}, 
with time reversal T. 
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In order to shed light on the peculiar ME effects 
in P orb ? we briefly show the group theory analysis 
and a derivation based on the Landau theory of phase 
transitions. [20| When working in the Cc space group with 
the symmetry operations {E, c}, the ferrimagnetic order 
leads to a lowered symmetry in the magnetic space group, 
with the presence of SOC. We define the order param- 
eter, M = — Sf + Sf, as a linear combination 
of Fe spins at A and B sites. Using the transformation 
rules given in Table HI we analyze, in the thermodynamic 
free energy, the possible ME coupling terms of the form 
P • M" 2 , which are invariant under symmetry operations. 



Fme = c aa P a Ml + c ah P a Ml + c ac P a M 2 c + c ca P c Ml 
+c ch P c Ml + c cc P c M 2 c + c aac P a M a M c 
+c cac P c M a M c + c bab P b M a M b + c bbc P b M b M c , (1) 

whereas the dielectric energy is traditionally written as: 
Fde = —P 2 /^Xi where the Qj, c^h coefficients and x 
(set as 1 in what follows) are constants. The minimum 
of F=F ME +F BE occurs when dF/dP a = dF/dP b = 
dF/dP c = 0, so that P is obtained. After assuming 
a simultaneous rotation of Fe spins in the be plane by 
an angle with respect to the c axis, setting M = 
M(0, sin#, cos (9), we derive 

PaiP) = ^(-c a6 + c ac )cos2#+^(c a6 + c ac ), 
P h (6) = Mi C66c sin2#, 

P c (0) = ^(-c c6 + c cc )cos2#+^(c c6 + c cc ). (2) 

Remarkably, the ^-dependence perfectly agrees with the 
experimentally observed property, P a oc sin 2 0.[d] Inci- 
dentally, we note that other ME systems often need two 
or more antiferromagnetic order parameters to cause ME 
effects [21], whereas here only one magnetic order param- 
eter M is relevant. 

Now, selecting the ^-dependent terms from Eq.(pQ), we 
obtain; 

Pme (0) = ^ (a P cos 2(9 + bp sin 2(9) , (3) 

where a P = (-c ab + c ac )P a + (-c cb + c cc )P c and 6 P = 
CbbcPb- Here, the ^-dependence of Fme is nothing but 
the MAE, i.e., we obtained an expression for the MAE 
as a function of P. This is consistent with what al- 
ready proposed by Rado et. al., where they assumed the 
electric field ^-dependence of the MAE coefficient 
Our additional advantage here is that we obtain the 0- 
dependence of P starting from the polar Cc symmetry, 
which was unknown at that time. Also note that none of 
the experimentally proposed centrosymmetric structures 
(e.g., P2/c or Prima) allows for TVf -induced P, because 
the inversion symmetry is not broken by the ferrimag- 
netic order (all Fe-B spins are parallel). Our findings re- 
lated to ME effects are therefore an indirect confirmation 
that, among those that were experimentally put forward, 
the Cc symmetry is the ground state of the CO-phase. 

DFT calculations — In order to quantitatively evalu- 
ate the P(0) behavior and to investigate the microscopic 
mechanism, DFT calculations were performed using the 
VAS P(22| code. Starting from a previous study by Jeng 
et al. [16], we used their optimized Cc structural param- 
eters and the GGA+£7[23[ approach, with U = 4.5 eV 
and J = 0.89 eV for Fe-d state (Fe-3(i 6 4 < s 2 electrons are 
treated as valence) . The 2x2x1 Monkhorst-Pack /c-point 
grid was used; other details are the same as in Ref. i 
A collinear ferrimagnetic configuration was set, with all 
Fe-B (Fe-A) sites as up-spin (down-spin) sites. By intro- 
ducing SOC self-consistently, we evaluated both i ) the 
MAE through the total-energy change and ii) the ME 
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FIG. 2: DFT results: (a) Change in on-site SOC energy 
(Psoc) and total energy (Ptotai) versus collinear spin direc- 
tion (0) in the be plane with respect to the b axis. The label 
xy/yz/zx denotes the each representative Fe 2+ -B site, (b) 
Change in P (/xC/m -2 ) versus with respect to P a =-3.86, 
Pb=0 and P c =4.89/xC/cm _2 at M//b. In order to compare 
with experimental results (Ref. |9|), the sign of AP a is re- 
versed and fitted to a function of a± sin 2 + a2sin 4 #, with 
ai = 1.24xl0~ 3 and a 2 =1.29xl0~ 5 (red solid line). 

effect through the change of P calculated by the Berry 
phase method, Q upon rotation of the M direction with 
respect to the crystalline axes. 

As for the MAE, our results are consistent with ex- 
perimental observations, i.e. of the magnetically hard 
a, intermediate b and easy c axes[9[ (with energy dif- 
ferences, E a — E\) = 4.2 meV/cell and — E c = 1.2 
meV/cell). Experimentally, the large afr-plane magnetic 
anisotropy is well known and used for the i^-field cooling 
technique, labeled "magnetic annealing" , which is gener- 
ally performed to grow large single-domain crystals below 
Ty.Q Following the procedure of Ref. [Tl, the "global" 
MAE can be decomposed into the "local" MAE, evalu- 
ated by the on-site SOC energy, Psoc [25] at each site. 
Figure [2] shows the "global" MAE and some selected "lo- 
cal" MAEs, upon rotation of Fe spins in the ab plane. 
Each local MAE shows a cos 20 curve, which can be ex- 
plained by SOC perturbation theory. On Fe 3+ (d 5 ) ions, 
a large MAE is not expected, due to the quenched or- 
bital states; on the other hand, on Fe 2+ (d 6 ) ions, we 
can apply the same discussion used for the Co 3+ (d 6 ) ion 
in BiCoOs[26j]. The perturbation theory predicts that, 
when the d^ y orbital is occupied, the single site MAE is 
given by E M ae — A 2 /2(l/A JT -4/A e J cos 20, where A is 
the SOC constant, Ajt is the JT splitting energy, A Gg is 



the energy splitting between (yz, zx)^ and states, and 
is the spin polar angle, giving energy minimum at 0=0 
(S//z). Our results therefore show that the local mag- 
netic easy axis is perpendicular to the occupied orbital 
plane (the latter being a hard plane); in other words, the 
local easy axis at xy/yz/zx site is the z/x/y axis, re- 
spectively. Indeed, in Fig. [2j the xy-\abe\ed MAE comes 
from the xy site, showing cos 20 curve with the minimum 
at M//z(~ c), whereas MAE from yz and zx site shows 
maximum nearby. The "global" MAE is the composition 
of these local MAEs, based on a delicate balance of (Sxy, 
Ayz, Azx) OO set; the subtle deviation of the local easy 
axis away from lattice vectors, due to the small tilting of 
the orbit als, also affects the final result. 

As for the focus of the paper, i.e. ME effects, we indeed 
obtained P ME (0) (cfr Figj2fb)), showing excellent agree- 
ment both with what expected from Landau theory in 
Eqj2] and with experimental results. Unfortunately, the 
experimental size of P ME has never been reported so far, 
so that our comparison cannot be quantitative. In our 
simulations, we got P ME < 20/iC/m -2 . This is roughly 
10 3 times smaller than P co ~ 5/iC/cm -2 and even one 
order of magnitude smaller than P ME =800/iC/m -2 in 
a representative ME system, TbMnOs [l2j]; however, it 
should be large enough to be measured, comparable with 
observed P ME =5/iC/ m - 2 in RbFe(Mo0 4 ) 2 .[27] Notably, 
the P5 trend put forward an exciting perspective: its sign 
can be switched by a 90°-rotation of the magnetic field 
(from -45° to 45° in Fig. [2(b)), pointing to a magnetic 
control of ferroelectric polarization. [28| 

In order to understand the mechanism driving ME ef- 
fects, let us examine how Fe-d orbital states are affected 
by the spin rotation in the be plane. Table HI1 summarizes 

TABLE II: 3d-orbital coefficients (in percentage, with spin 
states summed up) of the highest occupied d orbital and Fe-d 
occupancy rid on Fe 2+ B12 site with different SOC enhance- 
ment factors A (0=without SOC, xl=with standard SOC, 
and xl0=with the SOC term 10 times enhanced) for differ- 
ent M directions. 
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o '2 '2 '2 2 

xy yz zx Sz -r x -y 


n d 





98.85 0.02 0.37 0.76 


6.097 


xl M//b 
xl M//c 


97.97 0.4 0.45 0.78 0.41 
98.25 0.4 0.42 0.82 0.1 


6.083 
6.083 


xlO M//b 
xlO M//c 


51.74 17.25 16.82 1.94 12.24 
73.82 7.92 7.7 0.1 10.46 


6.108 
6.097 



the coefficients of a linear combination of <i-orbital states 
for the highest occupied d level [29| and the Fe-d occu- 
pancy at Fe-B12(xy) site when M//b and c. To highlight 
the effects of SOC, three calculations were done (keep- 
ing the ions fixed): (i) without SOC term (i.e. where 
the spin direction doesn't affect the orbitals), (ii) with 
standard SOC term, and (in) with SOC artificially en- 
hanced by a factor of 10. At first glance, d xy shows the 
largest component in all cases. In the extreme case with 
10-times-enhanced SOC, d yz and d zx orbitals have rela- 
tively large components, increasing to about 17% when 
M / /b, compared to only ~ 8% when M / / c. This result 
reflects the "inverse effect" of the perturbation theory 



discussed above, i.e. the orbital plane of the partly filled 
minority- spin state tends to be perpendicular to the im- 
posed spin direction, despite the fact that the structural 
JT distortion would favor a d xy state. The SOC also 
changes the occupancy of d electrons according to the 
spin direction. Although the change is negligibly small 
in standard SOC calculation, it is probably responsible 
for the slight change in P ME (0). Considering the above 
evidences, we propose the following scenario for magneto- 
electricity: i) SOC changes the shape of t^ g orbital when 
M is rotated (as shown in FigQJb)). ii) this in turn 
results in the anisotropic hybridization of the Fe-d state 
with surrounding O-p states, so that the gravity center of 
the d-electron (i.e. the Wannier center) is slightly shifted 
in one direction, hi) On the B12 site, the effect doesn't 
cancel out due to the lack of inversion-related ion in the 
unit cell, so that a purely electronic contribution appears 
in the net P orb . In other words, the fact that the underly- 
ing Cc symmetry lacks inversion is a necessary ingredient 
for ME effects to emerge. 

Conclusion. — We have shown that the origin of 
the ME effect in magnetite is "orbital-related" and P orb 
is induced by the on-site SOC term. This is clearly 
different from common conventional mechanisms (re- 
viewed in Ref . [3] ) . which invoke inter-site spin interac- 
tion, such as magneto-striction and DM interaction, in 
turn often connected to non-collinear spin-configurations. 
Rather, magnetoelectricity in magnetite emerge in a 
fully-collinear spin state with large magnetization, point- 
ing to an easier control of ferroelectric properties via mag- 
netic fields. We also note that the on-site P orb is not cou- 



pled with "structural" effects, such as piezo-electric or 
piezo-magnetic. Rather, the main ingredient which al- 
lows magnetite to be ME-active is basically the polar 00 
pattern of unfilled t2 g states. Indeed, no ionic displace- 
ments are involved in the rise of H-induced polarization 
and, in this sense, we can label this novel ME mecha- 
nism as purely driven by electronic (charge and orbital) 
degrees of freedom. Remarkably, our predicted trends 
show that the long-sought full control of polarization via 
a magnetic field can be achieved, with P5 switching sign 
upon a 90° rotation of H-field. We finally recall that 
the polar 00 is stabilized together with the polar CO 
in the Cc structure, through the CDW instability cou- 
pled with the Fermi-surface nesting below Tv-[30| These 
peculiar properties make magnetite rather unique; how- 
ever, we propose at least one candidate where the phe- 
nomenology might be similar: Ko.6FeF3,[3l| showing a 
non-centrosymmetric CO/OO pattern and where ME ef- 
fects should be explored, to confirm the existence of a 
novel class of OO-induced magnetoelectrics. 
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